if (!require("pacman")) install.packages("pacman")
pacman::p_load(dplyr, tidyr, stringr, purrr, fixest, ggplot2, broom)

# Load and filter data

dat <- readRDS("data/zij.rds") %>%
  filter(yob >= 1960 & yob <= 1975) %>%
  mutate_at(vars(one_of("DDR_kapitalismus", "vereinigung_dagegen", "DDR_reformsozial")), ~ . * 100)

# Create birth year groups for model

dat <- dat %>%
  mutate(birth_year_grp = ifelse(yob < 1971, "Pre1971", as.character(yob))) %>%
  mutate(birth_year_grp = factor(birth_year_grp, levels = c("Pre1971", "1971", "1972", "1973", "1974", "1975")))

# Define outcomes

outcomes <- c("DDR_kapitalismus", "vereinigung_dagegen", "DDR_reformsozial")

# Estimate models with birth year indicators

m2 <- feols(.[outcomes] ~ i(birth_year_grp, ref = "Pre1971"), data = dat, cluster = ~yob)
m2_women <- feols(.[outcomes] ~ i(birth_year_grp, ref = "Pre1971"), data = dat %>% filter(gender == "Women"), cluster = ~yob)
m2_men <- feols(.[outcomes] ~ i(birth_year_grp, ref = "Pre1971"), data = dat %>% filter(gender == "Men"), cluster = ~yob)

# Tidy and combine model results for plotting

tidy_feols_model <- function(model) {
  map_dfr(model, tidy, conf.int = TRUE, .id = "dv")
}

plot_data <- bind_rows(
  tidy_feols_model(m2) %>% mutate(sample = "All respondents"),
  tidy_feols_model(m2_women) %>% mutate(sample = "Women"),
  tidy_feols_model(m2_men) %>% mutate(sample = "Men")
) %>%
  filter(str_detect(term, "birth_year_grp::")) %>% # Keep only year coefficients
  mutate(
    birth_year = str_extract(term, "[0-9]{4}"), # Extract year from term
    birth_year = factor(birth_year, levels = c("1971", "1972", "1973", "1974", "1975")),
    dv = recode(dv,
      "lhs: DDR_kapitalismus" = "Support capitalism",
      "lhs: vereinigung_dagegen" = "Oppose reunification",
      "lhs: DDR_reformsozial" = "Support reformed socialism"
    )
  )

pd <- position_dodge(0.3)

# Generate plot

ggplot(plot_data, aes(x = factor(birth_year), y = estimate, group = sample)) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, color = sample), width = 0.0, position = pd) +
  geom_line(aes(color = sample), position = pd) +
  geom_point(aes(color = sample), position = pd, size = 2) +
  labs(
    x = "Year of birth (rel. to avg. of pre-1971 cohorts)",
    y = "Coefficient estimate\n(difference compared to avg. of pre-1971,\nin percentage points)",
    color = NULL
  ) +
  facet_wrap(~dv, scales = "free_y", ncol = 1) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    strip.text = element_text(size = 12),
    axis.title = element_text(size = 11)
  )
